Mechanistic role of alpha oscillations in a computational model of working memory

Brain oscillations are believed to be involved in the different operations necessary to manipulate information during working memory tasks. We propose a mechanistic role for the observed inhibition effect of the alpha rhythm based on its interference with the theta rhythm. Using the Lisman-Idiart model for multi-item working memory, we show that the interaction between these two oscillations is capable of creating a long lasting destructive interference that prevents the cyclic reactivation of neuronal ensembles and, as a consequence, memory maintenance. Additionally, to ensure robustness we propose a modular version of the model and implement oscillations as traveling waves. Using this model, we show that the interactions between theta and gamma determine the allocation of multiple memories in distinct modules, while the interference between theta and alpha disrupts the maintenance of the information already stored in them. The effect of alpha in erasing or blocking storage is robust and seems fairly independent of frequency, as long as it stays within the alpha range. This model helps us to understand why the alpha and theta oscillations, which have close frequency bands, could have opposite roles in working memory.


Introduction
Working memory is a putative memory system that incorporates many short-term information storage subsystems [1,2], and serves as an interface between perception, long term memory and action [3].In doing so, it contributes to higher cognitive functions such as reasoning, planning, decision-making and language comprehension [4].
Comprehensive models for a working memory system include specific components, or subsystems, to deal with different kinds of information and processes [2,5,6].In particular, for the short-term storage devices that compose the working memory system, most of the proposed models agree that some sort of persistent activity of the neurons involved in storing the information is necessary, but disagree on the underlying physiological mechanism.These models can be divided in roughly two categories: those that favor network mechanisms (e.g.instantaneous attractors due to short-term synaptic changes and continuous attractors or bump models) and models that favor single cell phenomena.
In the last category is the model proposed by Lisman and Idiart [7].According to this model, incoming information of a sequence of items to be memorized causes the ordered activation of item specific neural ensembles.The firing of these cells changes their intrinsic membrane properties producing a transitory excitability which peaks on a time scale of approximately 120ms after firing.In the absence of external inputs, these labeled cells can be reactivated by a nonspecific (non-informative) oscillatory input, provided that their period corresponds to the characteristic time of transient excitability.Therefore, they should be in the theta range.Another specific requirement of the model is feedback inhibition that prevents the synchronization of the different ensembles, so that the individual items are held active in time multiplexed fashion.A prediction of this model is a phase-amplitude coupling between the slow oscillations in the theta range representing the maintenance signal and the fast oscillations in the gamma range representing the firing of the stored memories.Subsequent studies have explored the theta-gamma mechanism in different models [8,9] (see [10] for a review).
Another oscillation that has been associated with an active role in cognitive processes, including working memory, is the alpha oscillation, which was the first rhythm observed in humans almost a century ago [11].Increases in alpha power were related to inhibited activity of areas during working memory maintenance tasks (see [12] for a review).During working memory scanning tasks, where subjects need to hold for a brief period of time several items in memory which are later tested, an increase in alpha power was observed during the retention period while a decrease was observed during the retrieval of information [13][14][15].An increase in alpha power related to the number of stored items was also observed [13][14][15].
A few theories have been used to describe possible roles for the alpha oscillation.The inhibition-timing hypothesis proposes that the peaks of alpha constrains time-windows of opportunities for the firing of neurons [16], with the increase of alpha power leading to a decrease of the windows size.The gating by inhibition hypothesis supposes that alpha can select the most relevant information to be processed through the blocking of irrelevant information routes [17].The oscillatory selection hypothesis suggests that the information selection could be accomplished by an entrainment between sensory stimulus and brain oscillations [18].In general, the main functional role addressed for alpha oscillations relates to the inhibition of irrelevant information.
Dippopa et.al. [19,20] proposed a computational model where brain oscillations act as functional operators of a working memory network, in other words, applying external oscillatory currents to the network allows it to store, maintain, prevent upload and erase information.In the model, oscillations in the beta-gamma range allow the upload of information, oscillations in the alpha range are responsible for erasing the content and preventing new content from being uploaded to the network, while oscillations in the theta range block upload but maintain the current content.However, in their current version of the model the issue of multiple item memory as well as the role of the oscillations phases are not fully explored.
Inspired by Dippopa et.al. [19], in this article we theoretically explore the possibility that oscillations in different frequency bands can play different roles in storage in a model with the characteristics of the Lisman-Idiart short-term memory model [7].More specifically, we want to explore whether increasing alpha power at the external oscillatory input can disrupt the multiple-item maintenance in the model therefore acting as a gating signal to erase and block memory storage.
The Lisman-Idiart model hinges on three essential elements: neural short-term excitability, concurrent external oscillation, and feedback inhibition.This paper, similar to the work of Dipoppa and Gutkin, [19], fits into the category of theory-based models outlined by Bassett et al. [21].In other words, we test our hypothesis exploring the potentiality of simple neural models, combining these elements without committing to specific biophysical mechanisms.
As for the network, we consider a modular version where the oscillations, responsible for the maintenance of the activity, sweep the network probing the modules as a traveling wave.
In the following sections we introduce the theta-alpha interference effect and the modular network model, discuss the results and present some conclusion.The mathematical descriptions of the neural models and the technical details about the computational experiments are left to the end in the Materials and Methods section.

The model
Here we present the modular version of the Lisman-Idiart theta-gamma model.The idea is to use space to improve the robustness of the model to noise as well as to enrich the possibilities of memory storage and manipulation.For simplicity we consider that the modules are distributed in a linear fashion (see Fig 1).In each module a short-term theta-gamma memory model is implemented as a network of excitatory neurons and inhibitory interneurons, that can perform multiple information storage through the cyclic reactivation of spatial firing patterns due to the combination of depolarizing currents after an action potential (ADP) and the excitatory drive by an unspecific oscillatory input on the theta-alpha frequency range .For most of the analysis, we utilize a current-based integrate-and-fire neural model (IF-Model) with an artificial resettable afterspike depolarization (ADP), as initially proposed by Lisman and Idiart [7].Additionally, as a proof of concept, we employ a more biophysically realistic model proposed by Rodriguez et al. [22] for L5 PFC pyramidal neurons, where the ADP is explained in terms of intrinsic calcium dynamics (high-threshold calcium currents and calcium-activated nonspecific cation current) in a Hodgkin-Huxley conductance based model (HH-Model).See Materials and methods.
In addition we consider that the oscillatory inputs are traveling waves in such a way that each module receives an oscillation with a slightly different phase.The evidence of travelling waves in the brain dates back almost to the first human EEG recording (see [23]), but its  importance has been recently highlighted by experimental work demonstrating that theta propagates in the hippocampus [24][25][26], and theta and alpha propagate in the neocortex [27].For a theoretical review see [28].
As in the original Lisman-Idiart Model, each module is capable of coding a number of overlapping firing patterns.But the modular version allows different items to be simultaneously stored in distinct modules that are scanned periodically by the travelling wave.It combines spatial and temporal properties such that the change in the oscillatory frequency input allows a much richer repertoire of operations for manipulating the information held in working memory.
Dipoppa et al [19] proposed four essential operations for a short-term memory model: load, maintain, block and erase.In this paper we show how each one of these operations can be performed in the context of our model, with the exception of the block operation that is not implemented independently but associated with the erase operation.
In other words, we consider that a network is blocked when it cannot hold information.The theta frequency will be responsible for allowing the load and maintain operations to take place and the alpha oscillations will perform the block/erase operations by interfering with the theta oscillations.

Results
Fig 2 illustrates the basic mechanism governing memory maintenance using both the IF-Model and the more intricate HH-Model.Both neuron models receive a sequence of three inputs: Memory maintenance.The model's memory behavior can be illustrated as follows: a neuron receives an informational input (square pulse) triggering firing, followed by a subthreshold network oscillation that, on its own, doesn't induce firing.If the amplitude of the network oscillation is too low the memory is not maintained.In (A), the neuron is an IF-model with stereotyped ADP, and in (C), it is the HH-model [22] with realistic calcium ADP and AHP.Panels (B) and (D) depict the impact of the oscillation phases due to the modularity of the network for both neuron models.https://doi.org/10.1371/journal.pone.0296217.g002an informational input sufficient to trigger firing, a subthreshold oscillation capable of sustaining firing in neurons that received the informational input but insufficient to cause firing on its own, and an oscillation unable to maintain firing.We observe that the more realistic neuron models, grounded in experimentally observed currents, replicates the behavior exhibited by the simpler model, Fig 2a.We employed the original set of parameters from Rodriguez et al. [22] for calcium currents (CaL, CAN) and calcium dynamics (Ca2), along with the parameters from Golomb and Amitai [29] for spike currents (Na, K), all without requiring further adjustments (refer to Materials and Methods).By supplying sufficient initial external current to activate after-depolarization (ADP), the HH-Model can be driven by an 8Hz oscillatory subthreshold current (see Fig 2c), demonstrating bistability through firing once per cycle during the oscillatory phase.The HH-Model proposed by Rodriguez et al. [22] exhibits an additional after-hyperpolarization (AHP) current, absent in our IF-Model.The significance of the AHP current lies in its ability to replicate spike irregularities observed in experimental data [22], this is less relevant in the context of the analysis conducted in this paper where the neurons are reactivated in a much lower frequency.Fig 2b and 2d demonstrate how firing can be temporally segregated in both models if the neurons are driven by oscillations with different phases.
Having established, at a qualitative level, the compatibility between the IF-Model and the HH-Model concerning the fundamental mechanism behind our memory model, we now restrict our analysis to the IF-Model.This choice is motivated by its fewer parameters, facilitating a clearer interpretation of the system.account for small asynchronies from the upstream inputs.The activities in each module are represented below with raster plots.The excitatory neurons are shown with a red gradient discriminating neurons that code different stimuli (A-D).Inhibitory neurons are shown in blue.On the bottom, the cumulative frequency for the firing of neurons representing each stimulus for each module are shown.The color code distinguishes each stimulus and the line style distinguishes the modules.Shaded areas represent the peaks of the traveling wave across the network.
In order to better understand the properties of the model after the introduction of the modular structure, we address the question of if there is an optimal delay between stimuli (or an optimal presentation frequency) for the correct allocation of items into different modules.We imagine four items (A through D) being presented to be held in memory at a rate f γ .This is translated in terms of a sequence of Gaussian input currents to the network (Fig 3 top).The Gaussian shapes of the inputs account for a certain degree of asynchrony in the upstream network representing the items and as a whole the input resembles a bout of gamma oscillations of frequency f γ (Fig 4a).The informational inputs are fed synchronously to all modules of the network, but they reach each of the modules on a different phase of the unspecific oscillatory input due its propagation speed (v osc ).The phase difference between consecutive modules is ψ osc = 2πf osc d/v osc where d is the physical distance between the center of two consecutive modules and f osc is the oscillatory frequency, and we make the simplifying assumption that all neurons in a given module are under the same phase.We define the input phase ϕ p,m of the p th item in the module m as the phase difference between the positive peak of the oscillation in that module and the item input time where ψ γ = 2πf oscd /f γ and ϕ i is the input phase of the first item in the first module (see Fig 4a).The input phase is positive if the stimulus anticipates the peak.Therefore, to evaluate the success of the load operation in different circumstances, three parameters are available: the presentation frequency of the items f γ , the phase difference between modules ψ osc and the input phase of the first item in the first module ϕ i .
We consider that the load operation is successful if items presented in a sequence are distributed in an orderly fashion among the different modules preserving sequence information.The easiest of possible orders is to have each item stored in a different module (Fig 3) so that the first item is represented in the first module, the second in the second, and so forth.The success of the load operation, therefore, depends on an optimal temporal match between the timings of the item inputs and the windows of opportunity defined by the positive phases of the oscillatory drive in each of the different modules.This involves adjusting the f γ , ψ osc and ϕ i parameters.We devised a criterion that optimizes f γ given a choice of ψ osc and ϕ i .for each set of (ψ osc , ϕ i ) and display the multi-item working memory parameter O s for the three initial cycles including the stimulation (left to right), showing a well defined region for a good working memory performance, exemplified by a rasterplot sample of three different conditions (blue, red and green).O s is a parameter defined to account for synchronization within modules and asynchronization between modules using only the statistical properties of the neuronal firing in order to evaluate the performance of the multi-item working memory storage (between 0 and 1).S1 The next issue is how well the network maintains information stored through time.In the Lisman-Idiart model for multi-item working memory, the combination of the oscillatory current in theta frequency and the membrane after-depolarization current creates a state of neuronal cyclic reactivation (Fig 5a), therefore neurons activated by a stimulus will be perpetually reactivated.
We consider that the maintenance operation is successfully accomplished if three requirements are met: i) all the neurons representing an item are firing with a good degree of synchrony during memory maintenance ii) each item is stored in a different module so neurons from different items fire asynchronously, iii) memory is reactivated in the first cycle of the unspecific oscillation.Using these criteria, we devised an order parameter (O s (z)) that measures the quality of the memory load and maintenance at a given cycle (z) (See Materials and methods).Higher values of O s (z) indicate better storage.Fig 5b shows that the four items correctly allocated are stably maintained through time by four different modules.
While the power of theta oscillations is correlated with the maintenance of working memory, the power of alpha is thought to actively inhibit irrelevant information on the same tasks.We propose that a possible mechanism for this inhibitory effect is the interference between existing oscillations in the alpha and theta range in the same network.We call it Oscillatory Interference Hypothesis.We consider that the combination of alpha and theta frequencies produces a beat, and the cyclic reactivation of the firing patters proposed in Lisman-Idiart ends up impaired and canceled, resulting in an erase operation on the working memory buffer (Fig 6a).In order to measure the correlation between the probability of alpha erasing the stored memories and three different parameters of the system (A α , α onset and f α ), we first calculate the average value of O s (z) parameter for the three cycles after alpha onset (z = 1 0 , 2 0 and 3 0 ), then binarized using the values of O s < 0.5 as memory successfully erased (P erase = 1) and O s > 0.5 memory not erased (P erase = 0) and plotted logistic regression fitting curves.Fig 5b shows that the inhibitory role of alpha was independent of the amplitude and the initial onset, but dependent on the specific frequency within the alpha band.

Discussion
We proposed an adaptation of the Lisman-Idiart model where the neurons are structured in spatial modules and the oscillations are traveling waves.The main reasons for these are (i) to increase the robustness of the model and (ii) to test the gating effect of different external oscillatory inputs.The original model segregated memories through fast-response inhibition, with sequential delays of the firing reactivations [7].As a result it was unstable to noise and prone to synchronization of neurons belonging to different items.It is important to note that underlying the concept of a multiplexing network, as proposed here, is Singer's Binding by Synchrony hypothesis [30].This hypothesis posits that the integral representation of an item is achieved through the synchronization of neural ensembles encoding its characteristics.For instance, the simultaneous representation in short-term memory of a blue circle and a red square is achieved through the specific synchronization of neurons corresponding to the accurate combination of shapes and colors.In simpler terms, blue neurons should synchronize with circle neurons, while red neurons should synchronize with square neurons.Consequently, synchronizing the ensembles of both items is damaging as it eliminates phase information, making it impossible to know the color associated with each shape.
With the introduction of network modularity and the oscillatory inputs as traveling waves the model becomes more robust.Recent evidence suggest that oscillations in the hippocampus and neocortex are indeed traveling waves [25][26][27].There is also overwhelming experimental evidence that the brain has a modular architecture, where neurons connect together forming micro local networks.The new version of the model, therefore, introduces not only phase but also space as part of the code.Since lists of items can come in any order we also assume that the modules have enough capacity to represent many items and that a specific item can be represented in different modules.This also solves the problem of storing similar items with overlapping neural representations since they will be active in different modules during multi-item storage.
In the simulated experiments we show that even if four stimuli are given equally to four modules, the excitatory dynamic created by a theta traveling wave allows the correct allocation of one memory item in each module.The optimal frequency to load a sequence of items is related to the speed of the traveling wave, that binds the theta and the gamma oscillations.This frequency is approximately the one in which the gamma period matches the difference in time between the phases of theta oscillations of two sequential modules being smaller as the global inhibition increases.In our model, the stimulus presentation should be interpreted as presented to the neural circuit and not the individual, being preprocessed by other brain regions and temporally compressed before being sent to the working memory system.Recent studies investigating theta and alpha oscillations as traveling waves indicate a phase difference of 2π rad/cm in the rat hippocampus [24], 0.1rad/cm in the human hippocampus [25] and 0.3 − 1.25rad/cm in the human neocortex [27], suggesting that our model would fit better on a human neocortical surface circuit, with a spatial scale of centimeters in and between modules.Our model indicates that there is a wide range of parameter allowing the theta and gamma oscillations to create the multi-item working memory storage.The blank area of the map of Fig 3b represents a more complex condition where more than one item is coded on the same module or the same item is coded by two different modules.Since every module has all the features of the Lisman-Idiart model, they could operate as local short-term buffers in case there is no interaction between modules.Regarding evidence about the role of the direction of propagation of the travelling waves in working memory [27], one possibility is that the coherence of the propagation of the travelling waves could recruit more or less modules, turning on the spatial features or maintaining just the local ones.
We use the same maintenance mechanism proposed on the Lisman-Idiart model, where stimulated neurons present afterdepolarization currents and the sum with the oscillatory theta input allows a cyclic reactivation.So, the same neurons once stimulated will fire again on each theta cycle.The same patterns activated by the stimuli A,B,C and D are repeated on each cycle for the modules M1, M2, M3 and M4.
The Oscillatory Interference Hypothesis showed to be a plausible mechanism for the blocking role of the alpha oscillations.Considering a condition where the alpha and theta oscillations compete for power (similar amplitudes).Alpha was able to disrupt the working memory performance in the model independently of its amplitude (Fig 6b , left panel).The effect also was relatively insensitive to the phase during the onset of alpha (Fig 6b, middle panel).As for the frequency values there is a sharp transition where, for a theta ' 8Hz, values of alpha frequency above 10Hz effectively erased the working memory buffer (Fig 6b, right panel).The measures taken used the combination of three cycles to ensure the time needed for the system to show long lasting and stable behaviour.
In our model alpha disrupts working memory by interfering with theta.The result is a beat profile that produces a long lasting period of low amplitude oscillations (Fig 7).Using a condition where both alpha and theta are synchronized at phase = 0 during alpha onset, it is possible to determine which values of alpha (and higher frequency bands) could produce a beat where the valley coincides with the peak of the ADP (I ADP > 0.85A ADP ). ) and theta frequency.This could explain the alpha frequency dependency found during our simulations.
Concerning the potential physiological validity of the derived properties within the IF-Model discussed here, we take initial steps towards a more realistic neural implementation.We showed that, at the fundamental level of memory maintenance, both the IF-Model and the more realistic HH-Model proposed by Rodriguez et al. [22] exhibit similar behavior.The intrinsic currents responsible for the HH-Model's after-spike depolarization (CaL, CAN) produce an effect akin to the IF-Model's more stereotyped ADP current.Both exhibit an asymmetric shape, characterized by a faster rise than decay, with amplitudes falling within the range of approximately 2.5 to 15mV.While Rodriguez et al. characterizes the ADP phenomenon as an after-burst depolarization, we consider every depolarization that emerges after each spike (See Fig 2).However, much exploration is needed within the broader parameter space of the HH-Model.

Conclusion
In this paper, we discussed a possible mechanism behind the effect that oscillations in the alpha range appear to have in cognitive tasks that demand that subjects disregard parts of the external stimuli.We do it in the light of the theta-gamma model proposed by Lisman-Idiart.In this model memory maintenance depends on two factors.The first is an intrinsic excitability caused by recent activity that tags neurons associated with a given memory.The second is an oscillatory input, in the theta range, that is sub-threshold for neurons that are not part of the memory but can drive the tagged neurons back to activity keeping them in oscillatory persistent activity, since there is a resonance effect between the time course of the excitability and the oscillatory input.According to the logic of the model, memory erasure could be accomplished by either eliminating the neural intrinsic response to firing or by disrupting the oscillatory input that refreshes the memory.We propose that an efficient way to disrupt the oscillation is by gating to the circuit an oscillation in the alpha range that will be superimposed over the existing theta oscillation causing amplitude modulations with the exact time scale to prevent the tagged neurons from refreshing their excitability.Although the same effect could be easily accomplished by reducing the power of the theta oscillations, in this paper we subscribe to the view that oscillations are natural attractors for the biological neural networks and preventing them may be more energetically costly than just combining them [31].In future work we aim to explore the interplay between oscillations and item similarity in long-term memory.This interaction has the potential to mitigate the impact of noise in the model.Noise disrupts firing patterns, causing neurons that should be firing together to desynchronize and synchronizing those that should not.Recurrent connections play a vital role in counteracting these tendencies, thereby preserving the integrity of neural ensembles that collectively represent the same memory.Once the parameter intervals for the basic IF-Model are established, a pivotal next step is to find equivalent parameter intervals in more realistic neural models, such as the one proposed by Rodriguez et al. [22], and subsequently compare them with physiologically realistic parameters.

Network
We consider a network composed by N neurons, of which N ex are principal excitatory and N inh are inhibitory interneurons.The network is divided in M spatial modules, with N a /M neurons of each type, a = ex, inh.The neural connectivity depends on the modules the neurons belong to, as well as the oscillatory inputs, with sequential oscillatory phase-differences producing the effect of a travelling wave.The stimulus input and the output occur simultaneously for all modules.

Connections
The network has the connectivity matrix presented in the S3 Fig, where the synaptic strengths are randomly generated by an uniform distribution between 0 and W type that depends on the kind of connection.
where W ab is the connection from the presynaptic neuron to the postsynaptic neuron and U(0, W type ) is a uniformly distributed random number between 0 and W type .W type can be Global (W EI , W IE ) or Modular (W EE , W EI , W IE ).

Integrate-and-fire neuron model
The neurons are modeled as current based integrate and fire, given by the equation where V is the membrane potential, V r is the resting potential, τ m is the membrane time constant and the last term is the sum over the i input currents.The potential is reset to a hyperpolarized value V reset after passing the firing threshold of V threshold , and stays unable to fire again for refractory time of t refractory .The total post-synaptic input received by the neuron "i" due to the firing of other neurons "j" in a given time "t" is where W ij is the synaptic weight, n j (t) is the number of spikes fired by the jth neuron up to time t, and t ðsÞ j are the spike times, and the individual input P is where H represent a Heaviside function.The principal excitatory neurons have an afterdepolarization potential that is reset for each new spike.
where t * i < t is the last spike of cell i before time t.The informational stimulus are modeled as Gaussian pulses where A inf is the input amplitude, t A the average time the stimulus was presented to the network and δ i 2 A indicates that only neurons linked to the information pattern A receive the inputs.The oscillatory inputs are given by the sinusoidal function where A osc is the amplitude, f osc is the frequency, ψ osc is a phase-shift creating a travelling wave effect and m i = 1, 2, . .., M is the module index of the ith neuron.The oscillatory power is modeled as all external oscillations come from the same source, meaning that the total power is a constrain condition for the system.

Noise
We consider only additive noise and introduce it in the simulations as a variability in the firing threshold for each neuron where is a normally distributed random variable with mean μ noise and standard deviation σ noise that is drawn independently for each neuron after each new spike.This is a strategy to decrease simulation time, valid for the IF-Model in the regime of low-frequency spiking (< 10Hz).Within this range, pure white noise added to the threshold is equivalent to more complete treatment of noise as long as its correlation decay time is significantly shorter than the theta oscillation period.

Neurophysiological realistic neuron model
Neuronal model adapted from L5 PFC pyramidal neurons from Rodriguez et al. [22].Neurons are Hodgkin-Huxley neurons described by The Leak current I L is written as The sodium current I Na is defined by The potassium current I K is defined by with The high-threshold calcium current I CaL is determined by the following equations The calcium-activated nonspecific cation current I CAN has the following equation with the activation x CAN depending on the intracellular calcium concentrations as The hyperpolarizing current I AHP is described by where de activation x AHP follows The calcium concentration dynamics is described by where The injection current is modeled in three modes as The parameters used are summarized in Table 2.

Metrics for loading performance
In order to evaluate the best frequency of stimuli presentation f γ , for loading information into the network given the parameters (ψ osc , ϕ i ), we simulated the loading cycle (zero cycle) of our network varying f γ between 100Hz and 33.33Hz.We then counted the number of activated cells for each item for each module n a,m where a 2 [A, B, C, D] and m 2 [M1, M2, M3, M4].We consider that, in a first approximation, a frequency is suitable for loading if the order of the stimuli is preserved in the modules.In other words, if the first item is the winner (the most active) in the first module, the second item is the winner in the second module, and so on.Mathematically a binary loading suitability can be written as where H[�] is the Heaviside function, the level g > 1 is a parameter that controls how bigger the winners must be (for instance, g = 2 indicates that the winner has to be at least twice the runner up), and the indexes i, j are numerical indexes representing items and modules considered here in equal number.We assume that the most suitable frequency (the best f γ ) for loading information for a given (ψ osc , ϕ i ) as the average When ∑ f ℓ(f|g) = 0 there is no suitable frequency, at level g, and the result is represented by a blank in Fig 4b.

Metrics for performance maintenance
We developed an order parameter that can account for the storage properties of the system.These are: neurons from a given stored item need to be synchronized, while stay asynchronized with neurons from other items.So, for the z th reactivation cycle and a stored item A, the measure for the synchronization within a memory item is where n A (z) is the number of active neurons in the ensemble that represents item A at the z th cycle, N A is the total number of neurons in the ensemble that represents item A, σ A (z) is the standard-deviation of the firing times of the n A (z) neurons, Δt is the time between the reactivation of two sequential items, β s is a control parameter that punishes the standard-deviation increase.
The measure for asynchronization between two memory items A and B, in the z th cycle is where

Metrics for erasing performance
In order to assess alpha's capacity of disrupting the stored memories, we performed 1200 simulations randomly varying three alpha parameters: α amplitude (A α , between 0.35 and 0.65 max), α onset (between 0 and 2π rad of θ ongoing oscillatory phase) and α frequency (f α , between 8 and 13Hz).To evaluate alpha's disrupting performance, we first calculate the average value of O s (z) parameter for the three cycles after alpha onset (z = 1 0 , 2 0 and 3 0 ), then binarized using the values of O s < 0.5 as memory successfully erased (P erase = 1) and O s > 0.5 memory not erased (P erase = 0) and plotted logistic regression fitting curves.

Analytical τ min curve
We computed the analytical τ min curve presented on Fig 6 as the half of the inverse of the beat frequency, which represents half of the beating period or the time to achieve the beating minimum.

Simulations and data analysis
We used Euler's method with step size dt = 0.01 ms for solving numerically the neurons differential equations.The simulations were written in C programming language and data analysis and graphic production were made with Python.

Fig 1 .
Fig 1. Modular working memory network.The modules are arranged linearly, each one composed by excitatory principal neurons (red circles) and inhibitory interneurons (blue circles).The informational inputs reach all modules synchronously while the unspecific oscillatory inputs (Osc 1 and Osc 2 ) are traveling waves and therefore, at a given time, the modules have oscillations with different phases.When in storage mode (maintain) a spatial firing pattern is reactivated for each module, producing the cyclic phenomena observed in the rasterplot of Fig 3. https://doi.org/10.1371/journal.pone.0296217.g001

Fig 2 .
Fig 2.Memory maintenance.The model's memory behavior can be illustrated as follows: a neuron receives an informational input (square pulse) triggering firing, followed by a subthreshold network oscillation that, on its own, doesn't induce firing.If the amplitude of the network oscillation is too low the memory is not maintained.In (A), the neuron is an IF-model with stereotyped ADP, and in (C), it is the HH-model[22] with realistic calcium ADP and AHP.Panels (B) and (D) depict the impact of the oscillation phases due to the modularity of the network for both neuron models.
Fig 3 illustrates the load operation, that allocates exactly one piece of information per module, when sequential stimuli are fed to the network.On the top we have the input currents of the four stimuli.The currents are Gaussian shaped to

Fig 3 .
Fig 3. Sequential allocation of items in different modules.On top: four Gaussian stimuli colored and labeled red for A, purple for B, blue for C and green for D. Middle: rasterplot of modules M1, M2, M3 and M4 during two cycles.The excitatory neurons are shown with a red gradient discriminating neurons that codes different stimuli (A-D).Inhibitory neurons are shown in blue.Shaded grey areas indicate the time windows of max neuronal excitability due to the theta wave of (>0.9 max).Bottom: cumulative frequency of the firing of neurons for each module (line style) and each stimulus (color).Label on the right, where a 2 [A,B,C,D].Parameters used: ψ osc = 0.9 rad/module, f γ = 50 Hz, f θ = 8 Hz, ϕ i = 0.8 rad.https://doi.org/10.1371/journal.pone.0296217.g003

Fig 4 .
Fig 4. Memory allocation on different modules.A) Definition of f γ , ψ osc and ϕ i .Four different stimulus are given equally to all four modules, with a given f γ frequency, modeled as independent Gaussian currents.ψ osc is the phase difference of the external oscillation of two sequential modules.ϕ i is the phase difference, in module M1, between the oscillations peak and the first input.B) Map of the f Best g for a given set of ϕ i and ψ osc when the oscillation is in the theta range f θ = 8Hz.Blank area represents the condition where the four stimulus were not correctly allocated to the four modules.C) f Best g dependency on ψ osc for three values of ϕ i .Shaded areas show the upper and lower values of f γ that loaded A,B,C and D in M1,M2,M3 and M4.D) ψ osc x ϕ i map of the working memory parameter O s using the f Best g during the first three cycles including stimulation (left to right).Three rasterplot samples are shown below in blue, red and green for three different conditions.The map is the mean of 50 repetitions.https://doi.org/10.1371/journal.pone.0296217.g004 Fig 4b shows a map where the color indicates the optimal values of f Best g , and where the blank area is the condition where the four stimuli are never allocated correctly (see Materials and methods).Fig 4c has the same information of Fig 4b, but represented as different curves of (f Best g � c osc ) for each ϕ i condition.In Fig 4d we take the f Best g Fig shows Fig 4c curve for all ϕ i conditions and S2 Fig shows Fig 4b and 4c using and alpha oscillation (12 Hz) instead of theta (8Hz).

Fig 5 .
Fig 5. Memory maintenance.A) Basic mechanism for maintenance at the neuronal level: the sum of the oscillatory theta input and the intrinsic afterdepolarization current is sufficient for a cyclic reactivation of the neurons once stimulated.B) Four items, one allocated to each module, are stably maintained in time for several theta cycles.The yaxis accounts for the fraction of active neurons of each item that during each cycle (n a (z)/N a , a 2 [A,B,C,D]) and the O s (z) parameter.https://doi.org/10.1371/journal.pone.0296217.g005 Fig 7 top colored curves shows two of the four beatings (red and green) with minimums (t 1 min and t 2 min coinciding with ADP's peak (ADP shown as the bottom curve, peak in bold) for theta = 8Hz.The analytical curve for the τ min is computed as half of the inverse of the beat frequency (see Materials and methods) and is showed as dashed gray line.The min and max alpha frequency that coincides with the ADP's peak is shown closed to the ADP's curve in black (10.4 Hz and 15.5 Hz, mean value 12.95 Hz).For other values of f θ , the inserted plot shows the dependence between the average value for alpha (f mean a

Fig 6 .
Fig 6.Memory erasing.A) Simulation's protocol: Load-Information is loaded in the network; Maintain-Information is maintained by reactivation for four oscillatory cycles; Erase-The onset of alpha oscillations occurs and storage measures are made on the next three cycles.B) Multi-item working memory performance measured by the O s parameter during the first three theta cycles (1',2' and 3') after the onset of alpha.The plots show dependency with the amplitude of alpha, the onset theta phase and the alpha frequency.Thick red lines shows the binned mean, where the markers are the bin's center and the shadowed area is the standard deviation.https://doi.org/10.1371/journal.pone.0296217.g006

Fig 7 .
Fig 7. Alpha inhibitory mechanism.A) Beat, in red, produced by an theta (8Hz) and alpha oscillation (10 Hz), dashed black lines, starting synchronized with phase = 0. B) Values of alpha (and higher frequency bands) could produce a beat where the valley coincides with the next two theta peaks, for theta = 8Hz.C) Mean possible alpha ( � f a ¼ f max a À f min a )for other values of f θ .https://doi.org/10.1371/journal.pone.0296217.g007

Table 1 .
Overview of parameters.Top: fixed parameters held constant throughout the simulations.Bottom: variable parameters uses more than one value or a range of values.